{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adaptive Derivative-Assembled Pseudo-Trotter ansatz Variational Quantum Eigensolver [1] (ADAPT-VQE)\n",
    "\n",
    "The Unitary Coupled Clusters truncated at the level of double excitations (UCCSD)\n",
    "$$ |\\psi^{UCCSD}\\rangle \\equiv e^{T(1)+T(2)}|\\psi^{HF}\\rangle, \\label{eq:1}\\tag{1}$$\n",
    "where\n",
    "$$ T(1)=\\sum_{ia} \\hat \\tau_{ia} = \\sum_{ia} \\tau_{ia} (a^{\\dagger}_{a} a_{i} - a^{\\dagger}_{i} a_{a}),$$\n",
    "$$T(2) = \\sum_{iajb} \\hat \\tau_{iajb} = \\sum_{iajb} \\tau_{iajb} (a^{\\dagger}_{a}  a^{\\dagger}_{b} a_{j} a_{i} - a^{\\dagger}_{i} a^{\\dagger}_{j} a_{b} a_{a}).$$\n",
    "is a natural ansatz for Variational Quantum Eigensolver since it is parametrized by the unitary operators and has a polynomial number of variational parameters. For the practical application on the quantum computer, it is desired to decompose unitary operator in Eq. (\\ref{eq:1}) into a product of unitary operators acting on a few qubits. This can be achieved via the first-order Trotter decomposition. Because of the variational flexibility of the ansatz (parameters are optimized after the Trotterization), one can sufficiently reproduce the UCCSD results performing only single Trotter step\n",
    "$$ |\\psi^{UCCSD}\\rangle \\approx \\prod_{ia}e^{\\hat \\tau_{ia}} \\prod_{iajb}e^{\\hat \\tau_{iajb}} |\\psi^{HF}\\rangle. \\label{eq:2}\\tag{2}$$\n",
    "The ADAPT-VQE algorithm suggests limiting the number of unitary operators in Eq. (\\ref{eq:2}) to lower the depth of the circuit and make the ansatz even more practical for quantum computing. This is achieved by the systematic extension of the ansatz by one excitation operator at a time, starting from the Hartree Fock state. Specifically, having specified the pool of all possible single and double excitation operators\n",
    "$\\{ \\hat A_1, \\hat A_2, \\dots, \\hat A_n \\}$ ($\\hat A_i \\in \\{ a^{\\dagger}_a a_i - a^{\\dagger}_i a_a; a^{\\dagger}_{a}  a^{\\dagger}_{b} a_{j} a_{i} - a^{\\dagger}_{i} a^{\\dagger}_{j} a_{b} a_{a} \\}$), at first one evaluates the gradient $$\\frac{\\partial E}{\\partial \\theta_m} = \n",
    "\\left.\\frac{\\partial}{\\partial \\theta_m} \\langle \\psi^{HF} |e^{-\\theta_1 \\hat A_1} e^{-\\theta_2 \\hat A_2} \\ldots\n",
    " \\hat H \\ldots e^{\\theta_2 \\hat A_2} e^{\\theta_1 \\hat A_1}| \\psi^{HF} \\rangle \\right|_{ \\forall \\theta_i=0} \n",
    " =\\langle \\psi^{HF} | [ \\hat H,\\hat A_m] | \\psi^{HF} \\rangle . \\label{eq:3}\\tag{3}$$\n",
    " Then ansatz is extended with the excitation operator,$\\hat A_k$, corresponding to the highest element of the gradient. Afterwards, all or only the new variational parameters are optimized. Continuing this process until the energy, the L$^2$ norm of the gradient or the maximum element is converged one recovers ADAPT-VQE [1]\n",
    "$$|\\psi^{ADAPT}\\rangle \\equiv (e^{\\hat \\tau_N}) (e^{\\hat \\tau_{N-1}}) \\ldots (e^{\\hat \\tau_2}) (e^{\\hat \\tau_1}) |\\psi^{HF}\\rangle, \\label{eq:4}\\tag{4}$$\n",
    "where $\\hat \\tau_p$ can be either $\\hat \\tau_{ia}$ or $\\hat \\tau_{iajb}$. It should be emphasized that the operator  pool is never drained and always contains all possible single and double excitation operators, thus allowing some excitation operator,$\\hat A_k$, enter the ansatz multiple times.\n",
    "\n",
    "[1] Harper R. Grimsley, Sophia E. Economou, Edwin Barnes & Nicholas J. Mayhall, *Natur. Commun.* 10, 3007 (2019)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Couldn't find cython int routine\n",
      "Couldn't find cython int routine\n"
     ]
    }
   ],
   "source": [
    "from qiskit import BasicAer\n",
    "#import warnings\n",
    "#warnings.filterwarnings('ignore')\n",
    "\n",
    "# setup qiskit_chemistry logging\n",
    "import logging\n",
    "from qiskit.chemistry import set_qiskit_chemistry_logging\n",
    "set_qiskit_chemistry_logging(logging.ERROR) # choose among DEBUG, INFO, WARNING, ERROR, CRITICAL and NOTSET\n",
    "\n",
    "# chemistry related modules\n",
    "from qiskit.chemistry import FermionicOperator\n",
    "from qiskit.chemistry.drivers import PySCFDriver, UnitsType\n",
    "\n",
    "#setup token for calculation on real device\n",
    "# from qiskit import IBMQ\n",
    "# provider = IBMQ.load_account()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Molecule Setup\n",
    "\n",
    "The hydrogen molecule (H$_2$) at the equilibrium bond length (0.735 Angstrom) is studied in this tutorial. In order to set up the molecular Hamiltonian one has to get the one- and two-electron integrals from a computational chemistry driver (PySCF in this case) employing the Hartree Fock calculation in STO-3G basis set. The information from the driver is saved in the ``molecule`` object. Then the Fermionic Hamiltonian is mapped to the qubit Hamiltonian using the ``parity`` mapping type as it further reduces the problem size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use PySCF, a classical computational chemistry software\n",
    "# package, to compute the one-body and two-body integrals in\n",
    "# molecular-orbital basis, necessary to form the Fermionic operator\n",
    "\n",
    "# chemistry related modules\n",
    "from qiskit.chemistry import FermionicOperator\n",
    "from qiskit.chemistry.drivers import PySCFDriver, UnitsType\n",
    "\n",
    "driver = PySCFDriver(atom='H .0 .0 .0; H .0 .0 0.735',\n",
    "                     unit=UnitsType.ANGSTROM,\n",
    "                     basis='sto3g')\n",
    "molecule = driver.run()\n",
    "num_particles = molecule.num_alpha + molecule.num_beta\n",
    "num_spin_orbitals = molecule.num_orbitals * 2\n",
    "\n",
    "from qiskit.aqua.operators import Z2Symmetries\n",
    "from qiskit.aqua.operators.legacy.op_converter import to_weighted_pauli_operator\n",
    "ferOp = FermionicOperator(h1=molecule.one_body_integrals, h2=molecule.two_body_integrals)\n",
    "\n",
    "# Build the qubit operator, which is the input to the VQE algorithm in Aqua\n",
    "map_type = 'PARITY'\n",
    "qubitOp = ferOp.mapping(map_type)\n",
    "# adapted this to work with master branch!\n",
    "qubitOp = Z2Symmetries.two_qubit_reduction(to_weighted_pauli_operator(qubitOp), num_particles)\n",
    "num_qubits = qubitOp.num_qubits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial State and Optimizer\n",
    "Set the initial state to the Hartree Fock state, $| \\psi^{HF} \\rangle$, and choose the optimizer for the VQE algorithm (in this case, limited memory BFGS from scipy)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup the initial state for the variational form\n",
    "from qiskit.chemistry.components.initial_states import HartreeFock\n",
    "init_state = HartreeFock(num_spin_orbitals, num_particles)\n",
    "\n",
    "# set the backend for the quantum computation=\n",
    "backend = BasicAer.get_backend('statevector_simulator')\n",
    "\n",
    "# setup a classical optimizer for VQE\n",
    "from qiskit.aqua.components.optimizers import L_BFGS_B\n",
    "optimizer = L_BFGS_B()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### UCCSD-VQE reference\n",
    "Perform the reference calculation employing the ``UCCSD`` variational form and ``VQE`` algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# setup the variational form for VQE\n",
    "from qiskit.chemistry.components.variational_forms import UCCSD\n",
    "var_form = UCCSD(num_spin_orbitals, num_particles, initial_state=init_state)\n",
    "from qiskit.aqua.algorithms import VQE\n",
    "algorithm = VQE(qubitOp, var_form, optimizer)\n",
    "result = algorithm.run(backend)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ADAPT-VQE run\n",
    "Now re-instantiate a new ``UCCSD`` object for the use with the ``VQEadapt`` algorithm, which is then internally made _adaptive_ by calling its ``manage_hopping_ops()`` method. In contrast to ``VQE`` one can additionally pass the ``threshold``and ``delta`` parameters to the ``VQEAdapt`` constructor. The  former defines the stopping criteria for the ansatz extension by setting the threshold for the maximum element of the gradient, the latter sets the accuracy of the gradient evaluation. In the current ADAPT-VQE implementation, the gradient is evaluated by the central difference approximation. One may also supply a custom pool of excitation operators to the optional `excitation_pool` keyword argument of the constructor. If it is omitted (like here) the pool is initialized to all possible single and double excitations in ``UCCSDAdapt``. Further keyword arguments are inherited from the ``VQE`` algorithm and are simply forwarded to any VQE calculations that are performed internally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "var_form_base = UCCSD(num_spin_orbitals, num_particles, initial_state=init_state)\n",
    "\n",
    "from qiskit.chemistry.algorithms import VQEAdapt\n",
    "algorithm = VQEAdapt(qubitOp, var_form_base, optimizer, threshold=0.00001, delta=0.1)\n",
    "result_adapt = algorithm.run(backend)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare results\n",
    "Print and compare the gates utilized in both parameterizations and the corresponding energies. Note the reduced number of `Evolution^1` gates in this simple example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------+---------------------+----------------------------------------------+\n",
      "|   Ansatz  |   Energy (Hartree)  |                    Gates                     |\n",
      "+-----------+---------------------+----------------------------------------------+\n",
      "|   UCCSD   | -1.8572750302023762 | OrderedDict([('Evolution^1', 3), ('u3', 1)]) |\n",
      "| ADAPT-VQE | -1.8572750302023793 | OrderedDict([('u3', 1), ('Evolution^1', 1)]) |\n",
      "+-----------+---------------------+----------------------------------------------+\n"
     ]
    }
   ],
   "source": [
    "from prettytable import PrettyTable\n",
    "table = PrettyTable()\n",
    "table.field_names = [\"Ansatz\",\"Energy (Hartree)\",\"Gates\"]\n",
    "cirq = var_form.construct_circuit(result.optimal_point)\n",
    "table.add_row(['UCCSD',str(result.eigenvalue.real),cirq.count_ops()])\n",
    "\n",
    "cirq_adapt = var_form_base.construct_circuit(result_adapt.optimal_point)\n",
    "table.add_row(['ADAPT-VQE', str(result_adapt.eigenvalue.real), cirq_adapt.count_ops()])\n",
    "print(table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also try running more elaborate examples like $H_3^+$ or $H_3^-$ which will lead to more obvious changes in the resulting gate operation counts. To do so, you can simply replace the string of atoms in the `PySCFDriver` instance further up top and set the charge to $\\pm1$. (e.g. `atom='H .0 .0 .0; H .0 .0 0.735; H .0 0.735 .0', charge=1` for a rectangular $H_3^+$)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
